Wavelet maxima curves of surface latent heat flux

ABSTRACT

The present invention introduces an innovative data mining technique to identify precursory signals associated with earthquakes. It involves a multistrategy approach that employs one-dimensional wavelet transformations to identify singularities in data, and analyzes the continuity of wavelet maxima in time and space to determine the singularities that could be precursory signals. Surface Latent Heat Flux (SLHF) data may be used. A single prominent SLHF anomaly may be found to be associated some days prior to a main earthquake event.

CROSS-REFERENCE TO RELATED APPLICATIONS

The present application is a continuation patent application that claimsthe benefit of nonprovisional utility patent application: Ser. No.11/108,115 to Cervone et al., filed on Apr. 18, 2005 now U.S. Pat. No.7,620,499, entitled: “Wavelet Maxima Curves of Surface Latent Heat Flux”and provisional patent application: Ser. No. 60/562,535 to Cervone etal., filed on Apr. 16, 2004, entitled “Wavelet Maxima Curves of SurfaceLatent Heat Flux Associated with the Greek Earthquake of 14 Aug. 2003,”both which are hereby incorporated by reference.

BACKGROUND OF THE INVENTION

Multisensor data available from airborne and spaceborne platforms havebecome widely used to study the changes of land, oceanic, atmosphericand ionospheric parameters, and their relation to various naturalhazards. For example, significant changes prior to earthquake eventshave been observed in Surface Latent Heat Flux (SLHF), Sea SurfaceTemperature (SST), water vapor and chlorophyll concentration. Suggestingthe presence of some interaction between the lithosphere and atmosphere,these observed changes have created an interest in using satellite-basedobservations to identify and study earthquake precursors.

Routine SLHF measurements can provide early warnings of an impendingearthquake. With respect to coastal earthquakes, anomalous SLHF peaksappear to consistently occur a few days prior to the main earthquakeevent. The magnitude of each peak tends to vary, while SLHF tends to behigher over oceans and lower over land. The origin of anomalous SLHFpeaks is likely to be related with the manifestations of the surfacetemperature in the epicentral region, which can be associated with thebuilding up of stress and movement along faults.

It is believed that temperature increases prior to an earthquake. Toshow changes in temperature, infrared (IR) wavelengths from a ModerateResolution Imaging Spectroradiometer (MODIS) sensor may be used. Theremay be various explanations, such as friction along a fault or fluidmovement, as to why the temperature may rise. In addition, SST may alsoincrease due to heat conduction. A rise in SST may cause oceanevaporation to increase, which in turn may raise anomalous SLHF peaksprior to a main earthquake.

Annually, SLHF contains a large number of maxima peaks, several of whichare more than 1 or 2 times above the standard deviation. These peaks canbe attributed to atmospheric phenomena, earthquakes or oceandisturbances. Problematically, it is difficult to identify the maximumSLHF peak as a precursor of an impending earthquake. The detection ofthe maximum SLHF peak is significant to alert and allow affected regionsto prepare for an impending earthquake. For example, had the maximumSLHF peak of the earthquake that struck the Indian Ocean in Dec. 26,2004 and caused the great tsunami disaster thereafter been detected,hundreds of thousands of people could have been forewarned and preparedfor evacuation. Hence, what is needed is a general methodology and modelto employ spatial and/or temporal analysis of wavelet maxima to identifysignals associated with earthquakes with precise continuity in time andspace.

Several commercial and research models have been developed for an earlywarning system that mainly uses past historical data, some of which canbe traced back to as far as the 5th century BC. For example, one modeluses past historical data that includes fracture zones calculated usinggravity fields, and changes in the electromagnetic field and tidalcycles to determine the occurrence of earthquakes. It further usesground based data. However, only one ground monitoring station isavailable. Thus, it precludes the possibility of real-time prediction.Therefore, what is also needed is a general methodology and device thatidentify signals associated with earthquakes in real-time.

BRIEF SUMMARY OF THE INVENTION

The present invention presents one aspect of identifying earthquakeprecursory signals by selecting a least one grid to generate a gridpath, analyzing the grid(s) in a time-series using wavelettransformation, identifying at least one local wavelet maximum,generating a sequence of local wavelet maxima, identifying at least onesingularity, measuring the geometrical space and time continuity of thesingularity along the grid path, and identifying an anomaly.

In yet a further aspect of the invention, at least one maximum line maybe propagated to reduce any potential numerical errors.

In yet a further aspect of the invention, a one-dimensional (1-D) realcontinuous wavelet transformation may be used to analyze a grid in atime-series.

In yet a further aspect of the invention, SLHF may be used.

In yet a further aspect of the invention, selecting a grid may be basedon the tectonics of a region, such as a two-dimensional (2-D) spacehaving dimensions width and height, a 2-D space having dimensionslatitude and longitude, evolutionary algorithms, heuristics, continentalboundaries, fault lines, and SLHF measurements.

In yet a further aspect of the invention, a grid path may be evaluatedby at least one criterion and tolerance. Each criterion used may bebased on at least one maximum length, at least one minimum spread, andat least one maximum anomaly.

In yet, a further aspect of the invention, geometrical measurements mayinvolve the introduction of a discontinuity in space and/or time forcompensating the rounding of errors in a wavelet transformation.

Additional objects, advantages and novel features of the invention willbe set forth in part in the description which follows, and in part willbecome apparent to those skilled in the art upon examination of thefollowing or may be learned by practice of the invention. The objectsand advantages of the invention may be realized and attained by means ofthe instrumentalities and combinations particularly pointed out in theappended claims.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

The accompanying drawings, which are incorporated in and form a part ofthe specification, illustrate an embodiment of the present inventionand, together with the description, serve to explain the principles ofthe invention.

FIG. 1 shows a method for identifying earthquake precursory signals asper one aspect of the invention.

FIG. 2 a shows an embodiment of an earthquake detector as per one aspectof the invention.

FIG. 2 b shows another embodiment of an earthquake detector as per oneaspect of the invention.

FIG. 3 shows extensional, compressional and strike-slip motions ofGreece's tectonics.

FIG. 4 shows a map of Greece regarding an Aug. 14, 2003 earthquake thatoccurred in Greece.

FIG. 5 shows wavelet analysis for SLHF data over Grid 28 for the year2003.

FIG. 6 shows wavelet analysis for SLHF data over Grid 29 for the year2003.

FIG. 7 shows wavelet analysis for SLHF data over Grid 30 for the year2003.

FIG. 8 shows wavelet analysis for SLHF data over Grid 39 for the year2003.

FIG. 9 shows wavelet analysis for SLHF data over Grid 31 for the year2003.

FIG. 10 shows wavelet analysis for SLHF data over Grid 22 for the year2003.

FIG. 11 shows wavelet analysis for SLHF data over Grid 23 for the year2003.

FIG. 12 shows wavelet maxima computed in wavelet transformations fordifferent grids, where the epicenter of the Aug. 14, 2003 earthquakeoccurred in Grid 31.

FIG. 13 shows results of possible paths of wavelet maxima computed inwavelet transformations.

FIG. 14 shows anomalies and time of occurrence for identified signals ateach grid location.

FIG. 15 shows yet another map of Greece, where the epicenter of the Aug.14, 2003 earthquake and the epicenter of the Mar. 1, 2004 earthquake aredenoted with stars.

FIG. 16 shows an example of wavelet analysis and corresponding maximacurves.

FIG. 17 shows wavelet analysis for SLHF data over Grid 37 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 18 shows wavelet analysis for SLHF data over Grid 38 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 19 shows wavelet analysis for SLHF data over Grid 39 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 20 shows wavelet analysis for SLHF data over Grid 40 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 21 shows wavelet analysis for SLHF data over Grid 49 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 22 shows wavelet analysis for SLHF data over Grid 41 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 23 shows wavelet analysis for SLHF data over Grid 32 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 24 shows wavelet analysis for SLHF data over Grid 33 from Mar. 30,2003 to Mar. 28, 2004.

FIG. 25 shows results of possible paths of wavelet maxima computed inwavelet transformations.

FIG. 26 shows the size of anomalies and their time of occurrence for thetwo precursory signals of the Aug. 14, 2003 earthquake (signal 1) andMar. 1, 2004 earthquake (signal 3).

DETAILED DESCRIPTION OF THE INVENTION

The present invention is a method, program storage device readable by amachine, and device that uses data mining techniques, including wavelettransformations and spatial and/or temporal continuity analysis ofwavelet maxima to identify earthquake precursory signals.

Specifically, the present invention may use wavelet transformations asdata mining tools by computing wavelet maxima that propagate from finerto coarser scales. These maxima may be used in identifying stronganomalies in data. Amongst these maxima, those showing continuity inboth space and time may be assumed to be potential precursors forearthquakes. Space continuity may be defined as: detected anomalies canbe distributed in space according to a precise geometry conforming to aregion's geological settings. Time continuity may be defined as:detected anomalies can occur at the same time or within a short delay ofeach other.

The present invention may be applied to SLHF data to study twoearthquakes that occurred in Greece: one on Aug. 14, 2003 and another onMar. 1, 2004. Significant SLHF anomalies may be found to be associatedto each earthquake event. While these two earthquakes were chosen merelyas examples since they occurred in the same region and within a shorttime span, it is important to note that the present invention is notlimited to a particular region or time frame for detecting earthquakes.

Referring to FIG. 1, earthquake precursory signals may be generallyidentified using a series of steps. These steps include: selecting atleast one grid S105, where a collection of grids can be used to generatea grid path; analyzing a time-series on the grid(s) using wavelettransformation S110; identifying at least local wavelet maximum for eachscale of the wavelet transformation S115; generating a sequence of localwavelet maxima S120; identifying at least one singularity using thesequence S125; measuring the geometrical space and time continuity ofthe singularity along the grid path S130; and identifying at least oneanomaly S135.

This method can be stored in a program storage device readable by amachine 205. Such device 205 can embody instructions executable by themachine to perform the method. Examples of program storage devicesinclude, but are not limited to, computers, compact discs (CDs), digitalvideo discs (DVDs), floppy disks, zip disks, hard drives, flash memorysticks/cards, random access memory (RAM), read only memory (ROM), etc.

To achieve this method, and other limitations below, an earthquakedetector 205 may be used. As illustrated in FIG. 2 a, the earthquakedetector 205 may have a grid selector 210, a wavelet transformationanalyzer 215, a local wavelet maximum identifier 220, a sequencegenerator 225, a singularity identifier 230, a geometric analyzer 235,and an anomaly identifier 240. It may also have a numerical errorsreducer 245, as shown in FIG. 2 b, to propagate at least one maximumline to reduce or ignore numerical errors. Examples of such deviceinclude, but are not limited to, personal digital assistants (PDAs), orits equivalent, mobile detector, centralized and/or localized detector,etc.

The grid selector 210 is capable of selecting at least one grid togenerate a grid path S105. The wavelet transformation analyzer 215 iscapable of analyzing wavelet transform in a time-series S110. The localwavelet maximum identifier 220 is capable of identifying at least onelocal wavelet maximum for each scale of the wavelet transformation S115.The sequence generator 225 is capable of generating a sequence of atleast one local wavelet maximum S120. The singularity identifier 230 iscapable of identifying at least one singularity using said sequenceS125. The geometrical analyzer 235 is capable of measuring space andtime continuity of at least one singularity along the grid path S130.The anomaly identifier 240 is capable of identifying at least oneanomaly S135.

Collectively, the method, program storage device and earthquake detectorcan share the following additional embodiments. First, any potentialnumerical errors, such as noise, may be reduced or ignored bypropagating at least one maximum line, or by considering only thelongest maxima lines. To help accomplish this feature, a numericalerrors reducer 245 may be used. Second, the grid(s) may be analyzed in atime-series using a 1-D real continuous wavelet transformation. Third,SLHF may be used. Fourth, to select a grid, it is preferable to use thetectonics of a region. Examples of tectonics include, but are notlimited to, a 2-D space having dimensions width and height (oralternatively longitude and latitude), at least one evolutionaryalgorithm, heuristics, at least one continental boundary, at least onefault line, and SLHF measurements. Fifth, the grid path may be evaluatedusing at least one criterion and at least one tolerance. The criterionmay be based on at least one maximum length, minimum spread and maximumanomaly. Sixth, the geometrical measurement may even include a step ofintroducing a discontinuity in space and/or time for compensatingrounding errors in a wavelet transformation.

The present invention may be applied to different types of spatialand/or temporal data. It is not meant to be bound to a particularspatial resolution or region in the world, or time sampling. Generally,earthquakes around the world are not randomly distributed. Rather, theytend to be concentrated in narrow, specific regions, such as along plateboundaries.

According to the United States Geological Survey (USGS), an earthquakeoccurs when there is a sudden slip on a fault that results in groundshaking and radiated seismic energy caused by the slip, or by volcanicof magmatic activity, or other sudden stress changes in the earth. Afault is a fracture along the Earth's crust in which the blocks of cruston either side have moved relative to one another parallel to thefracture. Fault fractures can be vertical or nearly vertical, orinclined. Where vertical or nearly vertical ones occur, blocks of crusttend to shift horizontally. Where inclined fractures occur, blocks ofcrust tend to shift vertically.

For purposes of exemplifying the usefulness and advantages of thepresent invention, the seismicity of Greece was considered.

Greece is located in the most seismically active region of theMediterranean and West Eurasian plate. This region is a part of thecollision zone between the Eurasian and African plates, whereearthquakes and volcanic eruptions are common. Two main different platetectonic boundaries southwest and east of Greece are the Hellenic Trenchand the Hellenic Arc. The Hellenic Trench is the largest area ofsubduction, in which the denser African plate goes under the less denseEurasian plate. In this region, magma rises from under the Earth's crustresulting in a large number of volcanoes. The Hellenic Arc is atransformation boundary, in which the African and Eurasian plates slideside by side, causing many fault lines which are responsible fornumerous earthquakes. However, there are no volcanoes in this region.

Typically, earthquakes are common in areas with high tectonicactivities. Greece has different types of tectonic regions, namelyextensional, compressional and strike-slip motions, as illustrated inFIG. 3. In these areas, earthquakes tend to exhibit large magnitude inthe compression and strike-slip zone, but small magnitude in theextension zone. Extensional motion, which can be coined as tensionalstress, is the pulling-apart type of motion. According to the USGS,tensional stress is the stress component perpendicular to a givensurface, such as a fault plane, that results from forces appliedperpendicularly to the surface of from remote forces transmitted througha surrounding rock. Compressional motion, which can coined ascompressional stress, is the squeezing type of motion. As the USGSdescribes, compressional stress is the stress component perpendicular toa given surface, such as a fault plane, that results from forces appliedperpendicularly to the surface or from remote forces transmitted througha surrounding rock. Strike-slip motion, which can be coined as sheerstress, is the stress component parallel to a given surface, such as afault plane, according to the USGS, that results from forces appliedparallel to the surface or from remote forces transmitted through asurrounding rock. According to the USGS, strike-slip fractures aretypically vertical blocks that have mostly moved horizontally.

According to the USGS, plate tectonics is a theory that combinesconcepts of continental drifting and sea-floor spreading. This theorysuggests that the Earth's rigid outer shell, known as the lithosphere,is broken into a plurality of oceanic and continental plates. Theseplates are in constant motion and have the ability to slide over theuppermost layer of the mantle, known as the asthenosphere. Currently,there exist seven major plates, which are subdivided into smallerplates. Each is approximately 80 kilometers thick; each is in constantmotion relative to one another. The rate of motion ranges between 10 to130 millimeters per year. The patterns of movement are neithersymmetrical nor simple.

Greece is known as a highly seismic and rapidly deforming region. OnAug. 14, 2003, an earthquake with a magnitude of 6.4 on the Richterscale occurred approximately 40 km off the Ionian island of Lefkada inWestern Greece. Fifty people were injured, and the region suffered heavydamages. The epicenter was approximately located at 39.18° N, 20.74° Ein the Ionian Sea with a focal depth of 10 km. Following this mainearthquake were two strong aftershocks, including a series of smalleraftershocks.

On Mar. 1, 2004, an earthquake of magnitude of 5.7 occurred about 7 kmnorthwest of the city of Kalamata. The area suffered light to moderatestructural damages. Fortunately, no casualties were directly associatedwith this seismic event. The epicenter was approximately located at37.23° N, 22.24° E in the SW Peloponnesos, with focal depth of 7 km. Theepicenter for this earthquake is located approximately 13 km north ofthe epicenter of the Aug. 14, 2003 earthquake.

Generally, in about every two years, an earthquake having a magnitude of5 or higher on the Richter scale occurs within 60 km of both theseepicenters. Both FIGS. 4 and 15 show a map of northern Greece. Starsdenote the location of an epicenter. Also shown are the location of theplate boundary, fault lines and their type (i.e., extensional,compressional and transformation).

As explained by the USGS, the Richter scale is a magnitude scaledeveloped in 1935 by Charles F. Richter of the California Institute ofTechnology. This scale serves as a mathematical device to measure andcompare earthquake sizes. A magnitude may be determined from thelogarithm of the amplitude of waves recorded by seismographs.Adjustments can be included to account for variations in distancesbetween various seismographs and the epicenter of earthquakes.

Magnitude may be expressed in whole numbers and decimal fractions. Forexample, an earthquake having a magnitude of 5.3 may be deemed as amoderate earthquake, whereas one having a magnitude of 6.3 may be deemedas a strong earthquake. Since the scale is based on logarithmic scale,each whole number increase represents a ten-fold increase in measuredamplitude. Equivalently, as an estimate of energy, each whole numberincrease corresponds to a release of approximately 31 times more energythan the amount associated with the preceding whole number. Instrumentsusing the Richter scale can be carefully calibrated with respect to eachother so that the magnitude may be computed from the record of anycalibrated seismograph.

Data used may include SLHF data from Jan. 1, 1998 to Oct. 31, 2003 forthe region bounded by latitudes 32° N to 44° N and longitudes 15° E to30° E. In addition, data may also include SLHF data from Jan. 1, 1998 toMar. 28, 2004 for the region bounded by latitudes 33° N to 45° N andlongitudes 14° E to 28° E. SLHF data may be obtained from the ScientificComputing Division of the National Center for Atmospheric Research(NCAR).

The data set may be represented by a Gaussian grid having 94 lines fromthe equator to a pole with a regular 1.8° longitudinal spacing andprojected into rectangular grid having dimensions of 2° latitude by 2°longitude. The global database of various meteorological parameters isoften maintained by National Centers for Environmental Protection(NCEP). This database may be generated by taking into considerationmeasured values at various worldwide stations and data retrieved fromsatellites. The fluxes used in operational weather forecast models mayincorporate in-situ observations through an assimilation process.However, the data source may be a frequent change of assimilationmethodology and model resolution, which has be solved by NCEP'swell-known re-analysis procedure that incorporates a whole archived dataset into a single, frozen data assimilation system.

Other data used may include those that identify and plot the location ofplate boundaries. It is well-known in the art that this data consists ofthe best fitting Euler vectors, closure fitting Euler vectors and theglobal model NUVEL-1. Such model may be used to geologically describecurrent plate motions between 12 assumed rigid plates. Additionally,data identifying and plotting slab contours may also be used. Such datamay include contours of the upper edge of subducting slabs calculatedusing relocated hypocentres. Furthermore, data used to identify and plotmajor fault lines may also be used. Moreover, data used to identify andplot the types of fault lines may be used. This data has been developedin the PLATES project by the Institute of Geophysics at the Universityof Texas.

The present invention's multi-strategy approach may employ a 1-D realcontinuous wavelet transformation to discover singularities in atime-series for a particular grid, and a geometrical analysis of thecontinuity in time and space of detected singularities across severalgrids adjacent to the epicenter chosen according to the tectonics of theregion, such as continental boundaries of fault lines, etc. The approachmay be tested using SLHF data. However, it is not constrained to aspecific type of data. For example, the approach may use PrecipitableWater data. The approach employing SLHF data may be implemented using asoftware called CQuake using, but not limited to, C, Java, Perl, Java Rcode, statistical package R, mathematical package octave developed atthe University of Wisconsin in combination with the WaveLab librarydeveloped at Stanford University, and General Mapping Tool (GMT)collection of functions developed at the University of Hawaii.

Wavelet Transformation

For any given real value function φ with zero average

∫_(−∞)^(∞)ϕ(t) 𝕕t = 0,let

$\begin{matrix}{{{Wf}\left( {u,s} \right)} = {\int{\frac{f(t)}{\sqrt{(s)}}{\phi\left( \frac{t - u}{s} \right)}{\mathbb{d}t}}}} & (1)\end{matrix}$be the real continuous wavelet transformation of a function ƒ. Since φhas zero mean, the previous integral measures the variation of ƒ in aneighborhood of time (position) u of size proportional to the so-calledscale factor s>0.

The set (u₀, s₀) may be defined to be a modulus maximum if |Wƒ(u, s₀)|is a local maximum, i.e., if

$\begin{matrix}{\frac{\partial{{Wf}\left( {u_{0},s_{0}} \right)}}{\partial u} = 0} & (2)\end{matrix}$and if Wƒ(u, s₀) is strictly increasing at the left of u₀ or strictlydecreasing at the right of u₀. In other words, the isolated local maxima(along the position coordinate u) of the wavelet transformation Wƒ(u, s)for each scale s>0 may be identified.

A connected curve γ in the scale-time place may be called a maxima lineif (u, s)εγ implies (u, s) is a modulus maximum.

Since the modulus maxima are isolated at each scale s₀, for each (u₀,s₀), there may be at most one curve γ that locally intersects the point(u₀, s₀). However, there may be none.

A singularity t₀ of the function ƒ may be defined as a point for whichthe derivative

$\frac{\mathbb{d}f}{\mathbb{d}t}$is not defined. Modulus maxima may carry a significant degree ofinformation about the position of singularities. More particularly, itmay be possible to prove that for every singularity, t₀ of ƒ, it may bepossible to find a sequence of modulus maxima (u_(i), s_(i)) such thats_(i)→0 and u_(i)→t₀.

It is preferred that the word “sequence” be used in lieu of the words“maxima line.” In general, maxima lines may break up and stop beforereaching the smallest scales. This break may present a problem in usingmaxima lines to identify singularities. However, there is a choice of φthat may be immune to this problem. If φ is a derivative of a Gaussian,then all modulus maxima may belong to maxima lines that propagate tofine scales. It is preferable to proceed with this choice.

It should be noted that even though sharp, isolated singularities of ƒmay be expected to be found at the abscissa of end points of significantmaxima lines, such finding may not always be the case. There can bemaxima lines that converge to perfectly smooth points of ƒ.

A rigorous detection of the singularities may also require thecomputation of the Lipschitz regularity coefficient of ƒ at the abscissaof the end point of the maxima lines. However, it is preferable to useonly the detection of significant maxima line as a possible indicationof the presence of sharp singularities.

As depicted in FIG. 16, a function with a non-smooth point at t=717 andthe corresponding set of maxima lines may be shown. The example showsthat some points converge for s→0 to 717, but there are some maximacurves converging to smooth points of the curve.

When dealing with signals with significantly large numerical errors, thenumber of maxima lines may also be very large. One way of lessening theimpact of high frequency numerical errors, such as noise, is to selectonly the longest or most important maxima lines. As a way to accomplishthis goal, a parameter called “propagation factor” (identified with ξ)may be introduced, where its values may be 0≦ξ≦1. This parameter may beused to determine the fraction of the total number of scales S that amaxima line γ must intersect to be considered as significant, i.e.,

$\begin{matrix}{\gamma = \left\{ \begin{matrix}{significant} & {{{if}\mspace{14mu}{length}\mspace{14mu}(\gamma)} \geq {\xi\; S}} \\{notsignificant} & {{{if}\mspace{14mu}{length}\mspace{14mu}(\gamma)} < {\xi\; S}}\end{matrix} \right.} & (3)\end{matrix}$

The length of γ may be computed from finer scales (high frequencies) tocoarser scales (low frequencies). It is preferable to consider assignificant the maxima lines that propagate at least up to scale≈ξS. Itshould be noted that S and ξS may be dependent on the specificdiscretization of the real continuous wavelet transformation. Whilethere may be an unavoidable degree of ambiguity in the choice ofdiscretization parameters, it is not likely to significantly affect thequality or robustness of the results.

While maxima lines were discussed in the context of continuous scalevariation and continuous position variation, in practice, theidentification of maxima lines may be performed in a discrete setting.Hence, there could be errors that prevent the very last point of themaxima line to converge to the probable singularity. Therefore, to havea more robust estimate of the time at which singularity occurs, thefollowing formula may be needed:

$\begin{matrix}{{g\left( {\gamma,n} \right)} = \frac{\sum\limits_{i = 1}^{n}\gamma_{i}}{n}} & (4)\end{matrix}$This formula tends to correspond to the average of the values at theabscissa of the last n end points of the maxima line, where n is aparameter either fixed or proportional to the number of scales, andγ_(i) is the abscissa value of γ at scale i. In cases having daily data,the calculated value g (γ, n) may correspond to the approximate day whenthe anomaly occurred. Experiments indicate that the average of the last10 values (n=10) can yield accurate determination of days correspondingto detected singularities. Yet still, it is important to note thatsmaller or larger values of n may be affected by the overall trend or bynumerical errors introduced by a signal's high frequency components.

The final result of the wavelet analysis is typically a characteristicvector of the same length as the original signal, where all values maybe set to zero except for those corresponding to the calculated valuesof g (γ, n) for all significant maxima lines γ. Such points may be setequal to the length of γ, or to the propagation of the maxima line. Theresulting characteristic vector may indicate both the time when asingularity was detected, and the singularity's weight as defined by thepropagation of the corresponding maxima line.

It should be noted that 1-D wavelet transformations are preferable inwavelet analysis over 2-D wavelet transformations because the latter maynot be able to detect the geometry of a signal or fault line. Moreover,similar 1-D construction along a continental boundary or fault line maystill be needed. Hence, the algorithm associated with 1-D measurementsmay be limited to avoid unnecessary complications.

Furthermore, 1-D wavelet transformations may effectively filter outpeaks due to small high frequency variations, which tend to beassociated with very short lines γ. Peaks caused by seasonal trends mayalso disregarded because the length of γ can be computed either fromfiner to coarser scales or higher to lower frequencies. Therefore,depending on propagation factor ξ, it may not be necessary if a line γpropagates all the way to the coarser scales or lowest frequencies,which can be affected by seasonal trends.

Continuity of Detected Singularities in Space and Time

The previously described wavelet analysis may be applied to atime-series of several grids corresponding to an epicentral region. Theresult of the analysis may identify several anomalous peaks, which maybe caused by earthquakes or by atmospheric or oceanic perturbations. Forexample, SLHF is an atmospheric parameter that is directly correlated tothe evaporation of water on the surface. SLHF may be affected by changesin land surface temperature and SST. The origin of an increase ofdecrease of surface temperature can be caused by atmosphericperturbations, such as strong winds, precipitation, intense cloud cover,or geological phenomena.

Such disturbances may be filtered starting with an assumption thatseismic activity is a large scale phenomena that manifests with aprecise geometry conforming to geological properties of a region. Peakscaused by earthquakes may be discriminated from peaks caused byatmospheric perturbations by using the concept of geometrical continuityunder space and time constraints. Not only may anomalous peaksassociated with earthquake events appear over a large area, but they mayalso follow a precise geometry as determined by a region's geologicalconditions. This characteristic generally means that a significantanomaly in one grid may be related to significant anomalies in othergrids if they all follow a precise geometrical path conforming togeological settings. Examples of such settings include continentalboundaries, fault lines, orography, sea depth and other characteristicsthat may help play a role on how anomalies might spread around anepicentral region. A collection of grids that satisfy these geometricalconstraints may be called “a grid path.”

It can be helpful to include a geometrical continuity of grid paths as aspatial constraint. This continuity can be used to single out any signalrelated to earthquakes. Furthermore, anomalies over a grid path may needto appear within a short period of time, such as a day or two. Thispresence may provide an additional constraint of the continuity of time.

Generally, spatial and temporal constraints can be used to “link,” in asuitable way, a singularity detected in each single grid to discriminatesignals most likely due to an impending earthquake from signals that arelikely due to other phenomena. As one aspect of creating this “link,”grid paths may need to be selected. Justification of the selection maybe necessary because a selection of contiguous grids defining the gridpath can greatly affect the results of the methodology employed.

However, the behavior of parameters involved may not respect fully anycontinuity along a suitably chosen grid path. This effect may be partlydue to atmospheric disturbances (such as wind, precipitation, etc.) ordiscontinuities in the sub-surface. Therefore, taking into account suchfactors, alternative grid paths should be selected.

Selecting grid paths may be viewed as a search in a 2-D space, where thedimensions are width and height. Alternatively, the dimensions can bebased upon longitude and latitude. Selection may be done using aspecialized search technique, such as evolutionary algorithms, orheuristics that can take into account different factors, which mayaffect the signal shape. Grids may even be selected by simply takinginto account continental boundaries or fault lines. Moreover, thetime-series measurements tend to be sufficiently smooth in space andtime. Such smoothness can potentially make the grid selection robustunder small perturbations of the grid path. The robustness factor maybecome more prevalent if higher resolution data is used.

Once the grid path is identified, the previously described wavelet-basedalgorithm may be applied and/or repeated for the time-series of selectedgrids to determine which signal is most likely associated with anearthquake.

The wavelet analysis may be performed using one of two ways: (1) astepwise format or (2) a smoothing format.

In the first way, results of a wavelet analysis for different grids maybe combined in a two dimensional matrix N(n×m), in which the row(s) n,correspond to time, and the column(s) m correspond to grid(s) of thegrid path at which the wavelet analysis may be performed. The values mayeither be the propagation length of significant maxima lines or zero ifno lines were detected in this particular point in space and/or time. Acontinuity line χ may be a signal that propagates in space and timeacross all non-zero cells of N for at least minLength columns and withinmaxDisc rows, where minLength and maxDisc are user defined parameters.

As one embodiment of generating continuity lines, a recursive algorithm,which take into account the following constraints, may be used:

-   -   1. C1: The minimum length of χ ought to be equal to or larger        than the parameter minLength.    -   2. C2: The time discontinuity can be at most ±maxDisc.    -   3. C3: The space discontinuity can be at most one consecutive        grid, but cannot be at either the first or last grid.

These constraints and derived formulas may be formalized using thefollowing:

$\begin{matrix}{{{\min\left( {\chi_{i}} \right)} \geq {\min\;{Length}}}{where}} & (5) \\{{\chi_{i}} = {{\sum\limits_{j = 1}^{m}A_{i,j}} - d_{i}}} & (6)\end{matrix}$is the length of a continuity line χ at time i, where

$\begin{matrix}{A_{i,j} = \left\{ {\begin{matrix}1 & {p_{i,j} > 0} \\0 & {p_{i,j} = 0}\end{matrix}\mspace{14mu}{and}} \right.} & (7) \\{p_{i,j} = \left\{ {\begin{matrix}\; & {{{if}\mspace{14mu} p_{{i + k},{j + 1}}} > 0} \\{N\left( {i,j} \right)} & \; \\\; & {{j - {\max\;{Disc}}} < k < {j + {\max\;{Disc}}}} \\0 & {otherwise}\end{matrix}\mspace{14mu}{and}} \right.} & (8) \\{d_{i} = {\overset{m - 1}{\sum\limits_{j = 2}}{B_{i,j}\mspace{14mu}{where}}}} & (9) \\{B_{i,j} = \left\{ \begin{matrix}\; & {p_{i,j} > 0} \\1 & {p_{{i + k},{j + 1}} = 0} \\\; & {p_{{i + k},{j + 2}} > 0} \\0 & {otherwise}\end{matrix} \right.} & (10)\end{matrix}$

The discontinuity in space and time may be introduced to compensate forrounding errors in the wavelet transformation. They may also helpcompensate for the approximation due to the averaging of the end pointsof significant maxima lines. Furthermore, they may help factor ingeophysical anomalies, such as wind, atmospheric perturbations and landorography.

Each grid path generated may be evaluated according to a multi-criterionevaluation function to determine the most important signals. Suchfunction may evaluate all signals using pairs of the form c, τ, where cis a criterion, and τ is a tolerance expressed in percentage. For eachcriterion c, it may compute the maximum value and select only the pathswithin tolerance of the maximum value. Typically, the paths whichsatisfy all criteria are considered to be significant. Three differentcriteria c may be used: 1) maximum length, 2) minimum spread and 3)maximum anomalies.

1. Maximum Length

This criterion evaluates the length of a path χ by summing the number ofgrids where a propagating signal may be present. As described earlier, adiscontinuity in a single grid may be allowed to compensate errors,which may be due to anomalies in fault line, geology or calculations.Taking into account such discontinuities, a parameter called“penalization coefficient” (PC) may be introduced to indicate a valuethat is subtracted for each discontinuity. For example, if thepenalization coefficient is equal to 0.5, a path of length 5 with nodiscontinuity is equivalent to a path of length 6 with twodiscontinuities. Alternatively, it is possible to add an increasinglylarger penalization coefficient to each new anomaly such that paths withnumerous discontinuities can be penalized more heavily.

2. Minimum Spread

This criterion evaluates the spread in time of a signal by consideringthe standard deviation for all paths and keeping only those paths withsmaller values. It may penalize signals that propagate across a faultline with a large discontinuity in time. For example, assume minDisc=2(days). Here, a signal would be moving two days forward in time across afault line. Hence, it should appear in grid 1 at time 100, grid 2 attime 102, grid 3 at time 104, etc.

3. Maximum Anomalies

This criterion takes into consideration any statistical significance ofanomalies associated with an earthquake event. Each wavelet maxima maycorrespond to a peak of an original signal. Its statistical significancemay be proportional to the size of the anomaly. Typically, largeanomalies are usually one to two sigmas above the mean value.

FIG. 12 shows the time of occurrence of wavelet maxima for differentgrids. It indicates points in space and time where significantsingularities in SLHF may be detected.

FIG. 13 shows the continuity of wavelet maxima in time and space.Significant signals are represented by horizontal lines. The day of theearthquake is shown with a dashed line. Peaks above one and two sigmasare respectively indicated with circles and triangles. The mostsignificant signal is shown with the thicker line.

Similarly, using the second way, results of a wavelet analysis fordifferent grids may be combined in a two dimensional matrix N (n×m), inwhich the row(s) n, correspond to grid(s) of the grid path at which thewavelet analysis may be performed, and the column(s) m correspond totime. The values may either be the propagation length of significantmaxima lines or zero if no lines were detected in this particular pointin space and/or time. A continuity line χ may be a signal thatpropagates in space and time across all non-zero cells of N for at leastminLength rows and within maxDisc columns, where minLength and maxDiscare user defined parameters.

As one embodiment of generating continuity lines, a recursive algorithm,which take into account the following constraints, may be used:

-   -   1. C1: The minimum length of χ ought to be equal to or larger        than the parameter minLength.    -   2. C2: The time discontinuity can be at most ±maxDisc.    -   3. C3: The space discontinuity can be at most one consecutive        grid, but cannot be at either the first or last grid.

These constraints and derived formulas may be formalized using thefollowing:min(|χ_(c)|)≧minLength  (11)specifying that the length of the continuity line χ at time c ought tobe greater than or equal to minLength, where

$\begin{matrix}{{\chi_{c}} = {\sum\limits_{r = 1}^{n}{\left( {A_{r,c} - {\sum\limits_{r = 1}^{n - 2}\left( {D_{r,c}*{PC}} \right)}} \right).}}} & (12)\end{matrix}$

This length may correspond to the number of detected anomalies A minusthe number of discontinuities for a chosen grid path D. Discontinuitiesmay be penalized by the parameter PC, which tends to indicate a valuethat is subtracted for each discontinuity. For example, if PC is equalto 0.5, a path of length 5 with no discontinuity is equivalent to a pathof 6 with two discontinuities. Alternatively, it is possible to add asan increasingly larger penalization coefficient to each new anomaly suchthat paths with many discontinuities can be penalized more heavily.

$\begin{matrix}{A_{r,c} = \left\{ \begin{matrix}\; & {{N\left( {r,k} \right)} > 0} \\1 & \; \\\; & {{c - {\max\;{Disc}}} < k < {c + {\max\;{Disc}}}} \\0 & {{otherwise}.}\end{matrix} \right.} & (13)\end{matrix}$

An anomaly A for grid r and time c may be detected if a wavelet maximumexists for grid r and time k. Time k may be defined as the interval ofthe original time c±maxDisc, and it can represent the maximum timediscontinuity allowed.

$\begin{matrix}{D_{r,c} = \left\{ \begin{matrix}\; & {N_{r,k} > 0} \\1 & {N_{{r + 1},k} = 0} \\\; & {N_{{r + 2},k} > 0} \\\; & {{c - {\max\;{Disc}}} < k < {c + {\max\;{Disc}}}} \\0 & {{otherwise}.}\end{matrix} \right.} & (14)\end{matrix}$

A discontinuity for grid r and time c may be detected if a waveletmaximum exists for grid r and grid r+2, but does not necessarily existfor grid r+1. The time k may be defined as given in Equation (14). Theremay be at most a discontinuity of a single grid.

Just as the stepwise format, the discontinuity in space and time mayalso be introduced to compensate for rounding errors in the wavelettransformation. They may also help compensate for the approximation dueto the averaging of the end points of significant maxima lines.Furthermore, they may help factor in geophysical anomalies, such aswind, atmospheric perturbations and land orography.

Each grid path generated may be evaluated according to a multi-criterionevaluation function to determine the most important signals. Suchfunction may evaluate all signals using pairs of the form Kτ, where K isa criterion, and τ is a tolerance expressed in percentage. For eachcriterion K, it may compute the maximum value and select only the pathswithin tolerance of the maximum value. Typically, the paths whichsatisfy all criteria are considered to be significant. Like the stepwiseformat, three different criteria K may be used: 1) maximum length, 2)minimum spread and 3) maximum anomalies.

1. Maximum Length

This criterion evaluates the length of a path χ by summing the number ofgrids where a propagating signal may be present. As described earlier, adiscontinuity in a single grid may be allowed to compensate errors,which may be due to anomalies along the selected grid path due torounding errors or inhomogeneity in the sub-surface.

2. Minimum Spread

This criterion evaluates the spread in time of a signal by consideringthe standard deviation (sigma) for all paths and keeping only thosepaths with smaller values. It may penalize signals that propagate acrossa fault line with a large discontinuity in time. For example, assumeminDisc=2 (days). Here, a signal would be moving two days forward intime across a fault line. Hence, it should appear in grid 1 at time 100,grid 2 at time 102, grid 3 at time 104, etc.

3. Maximum Anomalies

This criterion takes into consideration any statistical significance ofanomalies associated with earthquake events. Each wavelet maxima maycorrespond to a peak of an original signal. Its statistical significancemay be proportional to the size of the anomaly. Typically, largeanomalies are usually one to two sigmas above the mean value.

Finally, whether a stepping format or smoothing format is used, tocompare different earthquakes, the concept of anomaly A may beintroduced. Anomaly A may be defined as:

$\begin{matrix}{A = \frac{X_{i} - M}{STD}} & (15)\end{matrix}$where X_(i) is the SLHF at time i, M is the 30-day average for previousyears, and STD is the standard deviation of M.

Experiments

1. Stepping format

Referring to FIG. 4 again, using the Aug. 14, 2003 earthquake as anexample to determine the existence of a geometrical continuity ofanomalies over a continental boundary, grids 22, 23, 28, 29, 30, 31 and39 may be considered. The epicenter of the earthquake, shown with astar, is located in grid 31.

A wavelet analysis may be performed using Stanford University's Wavelablibrary and the University of Wisconsin's octave mathematical program.Table 1 shows wavelet parameters used for the analysis.

TABLE 1 Wavelet parameters of the numerical real wavelet transformation(RWT) used in experiment Function RWT Mother Wavelet Gauss Nvoice 20Scale  2 Octave 20 Propagation Factor ½

The wavelet analysis may be performed over the original signal for theyear 2003 for each grid. Data available from Jan. 1, 2003 to Oct. 31,2003 may be used. FIGS. 5-11 highlight results of wavelet analyses forSLHF data. These figures contain three parts. The first part shows atime-series for the original signal, a 30-day average for years1998-2002, and 1 and 2 standard deviations. The second part is agraphical representation of the previously described characteristicvector, which shows the time of significant wavelet maxima. The depth ofpropagation is shaded on a gray scale to emphasize those maxima thatpropagate to finer scales. The third part shows a graphicalrepresentation of the wavelet coefficients. Contour lines may be addedto increase readability and to stress the areas most likely to containmaxima lines. Table 2 shows information shared by FIGS. 5-11; Table 3identifies the location and grid used in each of these figures toconduct the wavelet analysis.

TABLE 2 Common SLHF Data Earthquake Date Aug. 14, 2003 Earthquake Day226 Wavelet Library WaveLab Function RWT (X, 20, 2, 10) PropagationFactor ½

TABLE 3 Location of Wavelet Analysis Figure Location Grid Number FIG. 539.05° N, 15° E 28 FIG. 6 39.05° N, 16.88° E 29 FIG. 7 39.05° N, 18.75°E 30 FIG. 8 40.95° N, 18.75° E 39 FIG. 9 39.05° N, 20.62° E 31 FIG. 1037.14° N, 20.62° E 22 FIG. 11 37.14° N, 22.5° E 23

A wavelet transformation for the epicentral region may be shown as inFIG. 9. It is possible to notice the sharp SLHF peak, which occurs abouttwo weeks prior to the earthquake, because the peak appears well abovetwo sigma (as represented by a dashed line).

Both FIGS. 12 and 13 show the analysis of spatial and temporalcontinuity. More specifically, FIG. 12 shows the local maxima of awavelet analysis where the X axis is the grid number and the Y axis istime expressed in days. FIG. 13 shows possible paths computed over theresults of the wavelet analysis. It is possible to notice that thereexists only one path that satisfies spatial and temporal constraints. Itis also possible to notice that it manifests about 12 days prior to theearthquake event. This signal may be considered as a strong precursor ofthe earthquake. The parameters used to compute the spatial and temporalcontinuity are given in Table 4.

TABLE 4 Continuity parameters used in the experiments Max Discontinuity20 Days Min Length 5 Grids Penalization Discontinuity 0.5 Tolerance 5%

FIG. 14 shows an SLHF anomaly for the only significant signal. Values onthe graph indicate the day when the anomaly was registered. These valuesare concentrated around day 214, corresponding to Aug. 2, 2003, just 12days prior to the earthquake event. The SLHF anomaly over the epicenterwas also found to be a maximum of 12 days prior to the main earthquakeevent.

In sum, the lessons learned from this example may be used to help detectearthquakes. The use of 1-D real continuous wavelet transformationcombined with the studies of spatial and temporal continuity of ananomaly over a continental boundary has discriminated a single signalamong the ones found. This signal occurred 12 days prior to theearthquake, and may be considered as a strong precursor. Such anomalousbehavior of SLHF may be directly associated with the occurrence ofcoastal earthquakes.

2. Smoothing Format

To exemplify how this format can use SLHF data to show a detectableanomaly prior to an earthquake, the Aug. 14, 2003 (having a magnitude of6.7) and Mar. 1, 2004 (having a magnitude of 5.7) earthquakes wereconsidered. Experiments may be carried out over the stronger earthquakeand validated using the same grid path over the weaker earthquake. Todetermine the existence of a geometrical continuity of anomalies over acontinental boundary, the grid path covering grids 32, 33, 37, 38, 39,40, 41 and 49 may be used. As shown in FIG. 15, the epicenter of eachearthquake, denoted with a star, is respectively shown in grids 41 and33. This figure also shows the location of the plate boundary, faultlines and their type (e.g., extensional, compressional, andtransformation).

Wavelet analysis may be performed using 365 SLHF data points, from Mar.30, 2003 to Mar. 28, 2004, for each of the grids in the grid path. Theresults of the wavelet analysis may be shown in FIGS. 17 to 24. Table 5shows the parameters used in the analysis.

TABLE 5 Wavelet parameters of the numerical RWT used in the experimentsFunction RWT Mother Wavelet Sombrero Nvoice 10 Scale 10 Octave  2Propagation Factor 12

Each of these figures may contain three parts:

1. The first part may show the time-series for the original signal, the30-days average for the previous five years 1998-2002, and the 1 and 2sigma for the 30-days average.

2. The second part may be a graphical representation of thecharacteristic vector as described in the previous wavelettransformation section, which shows the time when significant waveletmaxima may be detected. The gray scale generally indicates thepropagation depth to emphasize those maxima that propagate to finerscales.

3. The third part may include a graphical representation of the waveletcoefficients, and corresponding maxima lines. It is possible to noticehow the maxima lines converge to lines indicated in the previous part.

Table 2 shows information shared by FIGS. 17-24; Table 3 identifies thelocation and grid used in each of these figures to conduct the waveletanalysis.

TABLE 6 Common SLHF Data Earthquake Dates Aug. 14, 2003 and Mar. 1, 2004Data Analyzed Mar. 30, 2003-Mar. 28, 2004 Wavelet Library WaveLabFunction RWT (X, 10, Sombrero, 2, 10) Propagation Factor ½

TABLE 7 Location of Wavelet Analysis Figure Location Grid Number FIG. 1739.05° N, 13.12° E 37 FIG. 18 39.05° N, 15° E 38 FIG. 19 39.05° N,16.88° E 39 FIG. 20 39.05° N, 18.75° E 40 FIG. 21 40.95° N, 18.75° E 49FIG. 22 39.05° N, 20.62° E 41 FIG. 23 37.14° N, 20.62° E 32 FIG. 2437.14° N, 22.5° E 33

FIG. 22 shows the wavelet transformation for the epicentral region ofthe 2003 earthquake. A sharp SLHF peak well above two sigma (dashedline) is seen about two weeks prior to the earthquake. Similarly, FIG.24 shows the wavelet transformation for the epicentral region of the2004 earthquake. A sharp SLHF peak above one sigma (dotted line) is alsoseen about two weeks prior to the earthquake. The lesser intensity ofthis peak is likely due to the smaller magnitude of the 2004 earthquake.

Referring to FIG. 25, the analysis of spatial and temporal continuityillustrates significant paths generated using the results of the waveletanalysis, under specified time and space constraints. The x-axisrepresents time expressed in days; the y-axis represents the grid of thegrid path. The strength of the anomaly is shaded.

Out of a 365-day period, only three continuous paths may be found thatsatisfy space and time constraints. One of these paths occurred about 12days prior to the earthquake of Aug. 14, 2003. Another occurred about 14days prior to the Mar. 1, 2004 earthquake. Although the third signal maysatisfy the given constraints of space and time continuity, it appearsto be weaker than the former two. One reason is that the detectedsingularities for the third signal are not very strong. However, boththe stronger signals can be considered as strong precursors of the twoearthquakes. Table 6 provides the parameters used to compute spatial andtime continuity.

TABLE 6 Continuity parameters used in the experiments Max Discontinuity1 Day Min Length 8 Grids PC 2 Tolerance 20%

FIG. 26 shows SLHF anomalies for the two significant signals. The valueson the graph indicate the day when each of the anomalies was registered.It is possible to notice that in both cases, the largest anomalyoccurred either at the epicenter (earthquake of Aug. 14, 2003) or in thegrid immediately next to the epicenter (earthquake of Mar. 1, 2004).

In sum, the lessons learned from this example may be used to help detectearthquakes. The use of 1-D real continuous wavelet transformationcombined with the studies of spatial and temporal continuity of theanomaly over the continental boundary have shown SLHF anomalies asearthquake precursors. For both earthquakes, the observed precursorysignals are seen about two weeks prior to the main events. Theprecursory signals associated with the earthquake events may follow acontinuity in time and space, which can be used to discriminate othersignals due to different weather and atmospheric processes.

The foregoing descriptions of the preferred embodiments of the presentinvention have been presented for purposes of illustration anddescription. They are not intended to be exhaustive or to limit theinvention to the precise forms disclosed, and obviously manymodifications and variations are possible in light of the above teachingwithout departing from the scope of this invention and its broaderaspects. The illustrated embodiments were chosen and described in orderto best explain the principles of the invention and its practicalapplication to thereby enable others skilled in the art to best utilizethe invention in various embodiments and with various modifications asare suited to the particular use contemplated. For example, one skilledin the art will recognize that the present invention may be used in todetect seismic activity on Earth and planetary bodies, such as theEarth's moon.

A portion of the disclosure of this patent document contains materialwhich is subject to copyright protection. The copyright owner has noobjection to the facsimile reproduction by anyone of the patent documentor the patent disclosure, as it appears in the Patent and TrademarkOffice patent file or records, but otherwise reserves all copyrightrights whatsoever.

1. A program storage device readable by a machine, embodying a programof instructions executable by the machine to perform a process foridentifying earthquake precursory signals, the process comprising thesteps of: a. selecting at least one grid to generate a grid path; b.analyzing a time-series on the grid using wavelet transformation; c.identifying at least one local wavelet maximum for each scale of thewavelet transformation; d. generating a sequence of the wavelet maximum;e. identifying at least one singularity using the sequence; f. measuringgeometrical space and time continuity of the singularity along the gridpath; g. identifying at least one anomaly; and h. detecting anearthquake using the anomaly.
 2. The program storage device according toclaim 1, further including the step of propagating at least one maximumline to reduce numerical errors.
 3. The program storage device accordingto claim 1, wherein the process employs a one-dimensional realcontinuous wavelet transformation.
 4. The program storage deviceaccording to claim 1, wherein the grid is selected using the tectonicsof a region, including at least one of the following: a. atwo-dimensional space having dimensions width and height, b. atwo-dimensional space having dimensions longitude and latitude, c. atleast one evolutionary algorithm, d. heuristics, e. at least onecontinental boundary, f. at least one fault line, and g. Surface LatentHeat Flux measurements.
 5. The program storage device according to claim1, wherein the grid path is evaluated using: a. at least one criterion,including: i. at least one maximum length, ii. at least one minimumspread, and iii. at least one maximum anomaly, and b. at least onetolerance.
 6. The program storage device according to claim 1, whereinthe step of measuring geometrical space and time continuity furtherincludes the step of introducing a discontinuity in space and/or timefor compensating rounding errors in a wavelet transformation.
 7. Theprogram storage device according to claim 1, wherein the anomaly is aSurface Latent Heat Flux anomaly.